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We give a statistical-mechanical theory of stress transmission in disordered arrays of rigid grains 
with perfect friction. Starting from the equations of microscopic force and torque balance we derive 
the fundamental equations of stress equilibrium. We illustrate the validity of our approach by solving 
the stress distribution of a homogeneous and isotropic array. 
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Transmission of stress and statistics of force fluctua- 
tions in static granular arrays are fundamental, but unre- 
solved problems in physics . Despite several theoret 



ical attempts ]3j,|4| and a vast engineering literature 
the connectivity of granular media is still poorly under- 
stood at a fundamental level. In this Letter we propose a 
theory of stress transmission in disordered arrays of rigid 
cohesionless grains with perfect friction. A real granular 
aggregate (e.g. sand or soil) is a very complex object 
pL However, simple models are easier to comprehend, 
and extra complexities can always be incorporated sub- 
sequently. In our case the rigid grain paradigm provides 
a crucial starting point from which to appreciate the the- 
oretical physics of the problem. We model the granular 
material as an assembly of discrete rigid particles whose 
interactions with their neighbours arc localized at point- 
like contacts. Therefore the description of the network 
of intergranular contacts is essential for the understand- 
ing of force transmission in granular assemblies. Grain a 
exerts a force on grain (3 at a point lZ af3 = R Q + r a/3 . 
The contact is a point in a plane whose normal is n a P . 
The vector R a is defined by: 
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so that R Q is the centroid of contacts, and hence 
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where z is the number of contacts per grain and ^2 a 
means summation over the nearest neighbours. Hence 
R Q , r Q/3 and n are geometrical properties of the ag- 
gregate under consideration and the other shape spec- 
ifications do not enter. Friction is assumed to be infi- 
nite and the geometry is frozen after the deposition and 
can not be changed by applying or removing an exter- 
nal force on the boundaries. In a static array Newton's 
equations of integranular force and torque balance are 
satisfied. Balance of force around the grain a requires 
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where i — 1,2,3 are cartesian indices and g Q is the ex- 
ternal force acting on grain a. Further on in this Letter 
g is used also for the external forces at the boundaries. 

The equation of torque balance is 
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The centroid of the contact points need not coincide 
with the centroid of the forces e.g. the centre of mass of 
a solid grain, but we will assume it is so in order to keep 
the analysis simple so that we ensure that the macro- 
scopic stress tensor is symmetric, at least on average. It 
can be verified that, for the integranular forces in the 
static array to be determined by these equations, the co- 
ordination number z — 3 in 2-D and z = 4 in 3-D is 
required. In this paper we present the results for the 
2-D case only. The microscopic version of stress analy- 
sis is to determine all of the intergranular forces, given 
the applied force, torque loadings on each grain and ge- 
ometric specification of a granular array. The number of 
unknowns per grain is zd/2. Required force and torque 



equations give d 



d(d-l) 



constraints. The system of 



equations for the integranular forces is complete when 
the coordination number is z m = d + 1. Theory which 
confirms this observation has been proposed for periodic 
arrays of grains with perfect and zero friction |Q. It is 
clear that the coordination number z controls the connec- 
tivity of granular media. We will assume that z is indeed 
3 in 2-D, for this is surely the simplest situation, and one 
which is physically possible. The ultimate goal, however, 
is to determine the macroscopic stress tensor at every 
point of a granular array, given external loadings and ge- 
ometric specification. The macroscopic state of stress is 
a function of the distribution of contact forces. For any 
aggregate of discrete grains subjected to external loading, 
the transmission of stress from one point to another can 
only occur via the intergranular contacts. Therefore it is 
clear that the network of contacts determines the distri- 
bution of stresses within the granular array. The network 
of contacts is determined by the deposition history of the 
sample and the external loading on the boundaries. We 
define the tensorial force moment: 
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which is the microscopic analogue of the stress tensor. 
With C" = 0, will be symmetric. Our goal is to 
find a complete system of equations for the macroscopic 
stress tensor er^ , which is supported by the given network 
of contacts in the state of mechanical equilibrium. Given 
an assembly of discrete grains which is represented by a 
very complex network of contacts, we associate a continu- 
ous medium to have continuously distributed properties. 
Such spatial smoothing or coarse-graining can be accom- 
plished formally. To obtain the macroscopic stress tensor 
from the tensorial force moment to the macroscopic stress 
tensor we coarse-grain i.e. average it over an ensemble of 
configurations: 
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In the simplest cases of isotropic and homogeneous ar- 
rays this is not a problem. The difficulties appear when 
the array under consideration is anisotropic or inhomogc- 
neous. Within the confines of this paper we explore only 
the simplest cases. The number of equations required 
equals the number of independent components of a sym- 

At the 
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metric stress tensor tr^ — „j, u ,m ia „ 
same time, the number of equations available is d. These 
are vector equations of the stress equilibrium = gi 
which have their origin in Newton's second law. There- 
fore we have to find equations, which possess the 
information from Newton's third law, to complete and 
solve the system of equations which governs the trans- 
mission of stress in a granular array. Thus in 2-D there 
is one missing equation, and we derive it in terms of the 
geometry of the system. 

Given the set of equations (|3| — |5|) we can write the 
probability functional for the integranular force /- a ^ as 
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where the normalization, A/", which is a function of a con- 
figuration, is 



*~ 1= /n p {/"' 3 W Q/3 - 

Q,/3 



(9) 



where J T>f a P implies integration over all functions f a/3 , 
since all the constraints on / a @ have been experienced. 
We assume that the z = d + 1 condition means that the 
integral exists. 

We exponentiate the delta functions and thus intro- 



duce the set of conjugate fields Q™, 7", A" and 77, 
p{S9?} = jY[e iA Vf al3 V( a Vj a V\ a Vr 1 afi 
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where A is 
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The A a field term gives the symmetry of S°j. After inte- 
grating out the / and 77 Q/3 fields we find the following 
linear equation for the conjugate fields: 
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The idea of the conjugate fields method is to use these 
equations for the £ field in the stress probability func- 
tional, in order to derive the complete system of equa- 
tions for the stress tensor. The general solution of the 
above equation is a sum of the £° field which is the par- 
ticular solution and depends on 7, and C* which is the 
complimentary function 
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If we introduce the fabric tensor F," and its inverse 



M a : 
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13 
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we can rewrite the equation (13) in the following form: 





(16) 

which permits an expansion based on the first two terms, 



The probability of finding the tensorial force moment 
S?- on grain a is 

Pi^} = jj{s(s^-Y,f l a0 rf)p{fr P }'Dr P (10) 



i.e: 







km / j ' 
5 



7/) + 



(17) 



2 



The stress-force equation The next step is to in- where Fg is given by (na) and 



tegrate out the 7 field which gives us the stress-force 
equation: 

]T MgRfSS Y, M M R i° S f 3 = 9i ■ (18) 



So by expanding (3 quantities about a quantities we 
reach: 



After integrating out the £ Q * , 
the following equation for S£-\ 
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>" and ip a fields, we find 
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where K; 
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^ijki = \-fij n-j n k xtj ■') and gives a cor- 
rection to the standard equation of stress equilibrium 
Vj Uij = gi at the length-scale which is small compared 
to the size of the system. These corrections correspond 
to the presence of the second, third etc. nearest neigh- 
bours and topological correlations and must vanish in the 
k — > limit. 

The stress-geometry equation So far the well- 
known equations have been derived by using the informa- 
tion from Newton's second law. But we still have unused 
information from Newton's third law. By integrating out 
the £ * field we obtain the missing equations we are look- 
ing for. Let us consider that part of the eqn.(ll) which 
contains the £* field: 
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Note that although there are explicit forms generalis- 
ing ( p6[ ) in 3-D, these are more complex algebraically as a 
consequence of the higher coordination number. F£ and 
G-j will depend on configuration and averaging (|26|) is 
quite complex. The simplest array will have Q Q/3 or- 



thogonal to R Q , i.e. if R a/3 

qo/3 _ (y a P,-X a P 

can be written as 



(X a f 3 ,Y al3 ), then 
It follows, that F£ and G% 
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Then eqn.fPq) can be rewritten: 
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Counting the degrees of freedom in this equation we note 
that it can only give two (scalar) equations in 2-D and 



three in 3-D. Using R Q/3 and Q Q/: 



v a ^ r /3a 



can get these equations by projecting the vector equation 
equation into: 
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Thus if we are given S^ 2l the probability of finding 
Sf\ — S 22 is 
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Mathematically it is more convenient to introduce 
= (5a _ S^) 2 + (S^) 2 and determine the 
probability of finding Sf{ — given £ a . 
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It should be emphasised that the system under consid- 
eration is disordered and therefore Q 0/3 7^ (whereas 
for a honeycomb periodic array Q Q/3 = 0). Assuming 
Qj * gives rise to gradient terms we 
]) and (^||) by parametric variables 



as before that G? * 
can exponentiate 
4> a and ip a : 
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The mean values of 5*^ — S 22 and S^ 2 are zero, hence we 
predict, rather obviously, hydrostatic pressure. However, 
notice that we are able to predict the fluctuations away 
from hydrostatic pressure, and would do more on corre- 
lations if one could find a pathway to measuring them. 

Another approach to deal with the system of discrete 
equations (|]1| ^) is to solve it for S-j , and then average 
the solution. This way seems to be feasible at least for 
the simplest granular systems (e.g. isotropic or periodic 
arrays) and may provide deep insight into the origins of 
the non-gaussian statistics of stress fluctuations M. In 
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complex cases this can be accomplished in some approx- 
imation, or by using computer simulations. By applying 
Fourier transformation to (|l8|, ^6|) one can obtain £y (k) . 
The macroscopic stress tensor is obtained by averaging 
over the distribution of angles 8 a 
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array of fascinating theoretical and experimental prob- 
lems. Other important issues which are not addressed 
in this Letter are that of stress fluctuations and the re- 
sponse of a granular aggregate to external perturbations. 
In general, cohesionless granular materials are quasistatic 
or "fragile" Q, which means that they cannot support 
certain types of infinitesimal changes in stress without 
configurational rearrangements. 

In conclusion, our theory in its present form gives 
a simplified, but physical, picture of stress behaviour 
in cohesionless granular media. Further development is 
needed to make it a predictive tool which could be able 
to match experimental findings. 
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where |k| 2 = k\ + k\ and cr^-(r) = J a i j(k)e' 1 
By doing the inverse Fourier transformation one can see 
that the macroscopic stress tensor is diagonal. There 
must also be constraints on the permitted configurations 
(due to the absence of tensile forces) which are not so 
easily expressed, for they affect each grain in the form 
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which has not yet been put into continuum equations 
other than Deter > and Trcr > 0. 

DISCUSSION In this Letter we have derived the funda- 
mental equations of stress equilibrium: 
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In order to solve these equations one needs to know 
the geometric quantities K ijk i, Pi 3 k, T ikl and V i]k i- In 
practice details of the distribution of intergranular con- 
tacts are not known in advance, but should be obtained 
from the deposition history of the system or experimental 
measurements of two-body correlation functions. 

If the system is strongly anisotropic (i.e. there exists 
a preferred direction characterised by some ang le <fi) and 
tan 9 a has an average value tan <j>, then eqn.(|26|) becomes 
in the mean-field approximation 
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where <f> is the angle of repose. It is known as the Fixed 
Principal Axes equation |]3j, and has been used with no- 
table effect to solve the problem of the stress distribution 
in sandpiles. Explicit mathematical expressions for the 3- 
D case are more complex, and will be reported elsewhere. 
The issue of whether the derived system of equations ( |35| - 
|36| ) is robust against the inclusion of real friciton, soft- 
ness of grains etc. illuminates the existence of a whole 
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